métabo

DUBii 2021

Chargement des librairies

Attention select de mass (dépendance de mixOmics) prend le pas sur dplyr

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.1     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library("mixOmics")

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: lattice

## 
## Loaded mixOmics 6.14.1
## Thank you for using mixOmics!
## Tutorials: http://mixomics.org
## Bookdown vignette: https://mixomicsteam.github.io/Bookdown
## Questions, issues: Follow the prompts at http://mixomics.org/contact-us
## Cite us:  citation('mixOmics')

## 
## Attaching package: 'mixOmics'

## The following object is masked from 'package:purrr':
## 
##     map

library(FactoMineR)
library("factoextra")

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(skimr)

library(knitr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

Chargement des données brutes

chemin <- "../Raw_brut/"

#Echantillons
Sample <- read.csv2(paste0(chemin,"MI_Covariates.csv"))
Corres_var <- read.csv2(paste0(chemin,"MI_Covariates_corresp.csv"))
Corres_ID <- read.csv2(paste0(chemin,"Mb_Info_seq.csv"))

#Métabolites
NMR_CPMG <- read.csv2(paste0(chemin,"NMR_CPMG.csv"))
NMR_NOESY <- read.csv2(paste0(chemin,"NMR_NOESY1D.csv"))
List_cplte_CPMG <- read.csv2(paste0(chemin,"NMR_BUCKET_CPMG_LIST_assigned.csv"))
List_cplte_NOESY <- read.csv2(paste0(chemin,"NMR_BUCKET_NOESY_LIST_assigned.csv"))
List_CPMG <- read.csv2(paste0(chemin,"List_NMR.csv"))
List_NOESY <- read.csv2(paste0(chemin,"List_NOESY.csv"))


## Homogénéisation des noms de variables
##Attention package mass possède un select qui annule celui de dplyr
Sample_ID <- left_join(Sample,Corres_ID,by=c("DonorID"="DonorId")) %>% 
  relocate("Echantillon", .before = "DonorID")  %>% 
  dplyr::select(-"SeqDepth",-"DonorID")
rownames(Sample_ID) = Sample_ID$Echantillon

#Met la variable AGE en catégorie non hierarchisé et transformation en factor
Sample_ID <- Sample_ID %>% 
  dplyr::select(-"Echantillon") %>% 
  mutate(age = factor(AGE, levels = c(1, 2, 3, 4, 5), labels = c("[20-29]", "[30-39]", "[40-49]", "[50-59]", "[60-69]"))) %>% 
  mutate(SEX = factor(SEX)) 
Sample_ID$tabac <- factor(Sample_ID$tabac)



CPMG <- left_join(NMR_CPMG,Corres_ID,by=c("DonnorID"="DonorId")) %>% 
  relocate("Echantillon", .before = "DonnorID") %>% 
  dplyr::select(-c("DonnorID","Sample_type","VagueManip","Date","Spectre2","Exp","SeqDepth"))
rownames(CPMG) = CPMG$Echantillon
CPMG <- CPMG %>% 
  dplyr::select(-Echantillon)

NOESY <- left_join(NMR_NOESY,Corres_ID,by=c("DonnorID"="DonorId")) %>% 
  relocate("Echantillon", .before = "DonnorID") %>% 
  dplyr::select(-c("DonnorID","Spectre2","Exp","Sample_type","VagueManip","CV_Grps","CV_Grps6","CV_Grps7","CV_Grps8","CV_Grps9","SeqDepth"))
rownames(NOESY) = NOESY$Echantillon
NOESY <- NOESY %>% 
  dplyr::select(-Echantillon)

# Passage en matrice
Noesy <- as.matrix(NOESY)
Cpmg <- as.matrix(CPMG)




# ?mixOmics
# #selection des métabolites d'interêts. Voir pour transposer et ne selectionner que certaines col avec métabo nommés 
# List_CPMG1 <- t(List_CPMG)
# colnames(List_CPMG1) = List_CPMG1[1,]
# List_CPMG2 <-  t(List_CPMG1[3,])
# List_CPMG2 <- t(List_CPMG2)

Représentation des données CPMG

Histogramme

tCPMG_vector <- unlist(as.vector(t(CPMG)))
CPMG_vector <- unlist(as.vector(CPMG))

hist(CPMG_vector, breaks = 100)

hist(log10(CPMG_vector + 1e-8), breaks = 100)

df_row_sum <- apply(CPMG,MARGIN = 1, FUN= mean)
df_col_sum <- apply(CPMG,MARGIN = 2, FUN= mean)
#view(df_sum)
meanr <- mean(df_row_sum)
meanc <- mean(df_col_sum)
minr <- min(df_row_sum)
maxr <- max(df_row_sum)
minc <- min(df_col_sum)
maxc <- max(df_col_sum)
medr <- median(df_row_sum)
medc <- median(df_col_sum)


df <- data.frame(
  nb_ech = c(846,202),
  moy = c(meanr,meanc),
  median = c(medr,medc),
  min = c(minr,minc),
  max = c(maxr, maxc)
)
rownames(df) <- c("Echantillon","Bucket")
kable(df, align = "c", caption = "Tableau de moyenne des échantillons et buckets") %>%
  kable_styling(
    #bootstrap_options = "striped",
    full_width = F,
    position = "center")

Tableau de moyenne des échantillons et buckets

	nb_ech	moy	median	min	max
Echantillon	846	0.0059508	0.0059521	0.0058323	0.0060026
Bucket	202	0.0059508	0.0015731	0.0000588	0.1327695

#min(CPMG)
#max(CPMG)

BoxPlot

knitr::opts_chunk$set(fig.height = 3.5, fig.width = 8.5)

log10_tCPMG <- log10(t(CPMG) + 1e-8)
log10_CPMG <- log10(CPMG + 1e-8)

# par_ori <- par(no.readonly = TRUE) ## Store the original parameters in order to restore them afterwards
# par(mar = c(4, 6, 5, 1)) ## Adapt the margins of the boxplot in order to have enough space for sample labels
# par(mfrow = c(1,2)) ## Two side-by-side panes


boxplot(t(CPMG), 
        horizontal = TRUE, 
        las = 1, 
        main = "original values", 
        xlab = "value")

boxplot(log10_tCPMG, 
        horizontal = TRUE, 
        las = 1, 
        main = "log10-transformed", 
        xlab = "log10(value)")

Il y a trop de buckets pour que ce soit exploitable. Sachant que seuls 20% sont assignés à des métabolites, il faudrait peut être voir pour ne prendre que ceux là ??

Au vue des boxplot qui sont difficilement exploitables, je suis resté sur la transformation log10(CPMG + 1e-6) Avant transformation, il y avait en sorte 2 groupes : inf et supérieure 0.10 Après transformation cela semble moins marqué mais il reste une “bande” > -1 A mon avis quelque soit la transformation, cela restera. Dois-je continuer pour les centrer ??

Visualisation des ACP

Les 2 premiers axes expliquent plus de 62% de la variance exprimée

res.pca <- PCA(log10_CPMG, graph = FALSE)
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 60))

• En fonction des groupes :

Quelque soit le groupe étudié (sexe, tabac, âge), sur les axes 1 et 2, les différents points se superposent, il n’est donc pas possible de distinguer un groupe d’un autre. Visualisation montrer en fonction du groupe sexe Il y a un point complétement excentré (18%axe1 et 80% axe2) Homme 60-69ans exfumeur. Doit on garder cet individu (L31R124) ? C’est la plus petite moy (583e-5 et les autres >590e-5) max et min pas si différents des autres

fviz_pca_ind(res.pca, center=TRUE ,geom.ind = "point", habillage = Sample_ID$SEX, addEllipses = TRUE)

fviz_pca_ind(res.pca, geom.ind = "point", habillage = Sample_ID$SEX, addEllipses = TRUE)

fviz_pca_ind(res.pca, geom.ind = "point", habillage = Sample_ID$age, addEllipses = TRUE)

fviz_pca_ind(res.pca, geom.ind = "point", habillage = Sample_ID$tabac, addEllipses = TRUE)

Sur les axes 3 et 4 il est plus facile de distinguer les 2 sexes même s’il y a du chevauchement. Pour les 2 autres il n’est pas posssible de les discrimer. Il y a également 2 points completement excentrés à -20 et -50 axes 3 et 4. On les garde ? Pas discordant dans les autres métadatas

fviz_pca_ind(res.pca, axes= c(3,4), geom.ind = "point", habillage = Sample_ID$SEX)

# fviz_pca_ind(res.pca,  axes= c(3,4), geom.ind = "point",
#              habillage = Sample_ID$age, addEllipses = TRUE)
#
# fviz_pca_ind(res.pca, axes= c(3,4), geom.ind = "point", habillage = Sample_ID$tabac)

Cercle des correlations

fviz_pca_var(res.pca, col.var = "black", geom = "point")

Biplot des 40 buckets qui contribuent le plus à la variance expliquée

La contribution de ces 40 buckets est entre 0.7 et 0.75

fviz_pca_biplot(res.pca, repel = TRUE, axes = c(1,2),
                select.var = list(contrib = 40),
                                  geom.ind = "point",
                                  fill.ind = Sample_ID$SEX,
                                  pointshape = 21, pointsize = 2,
                                  palette = "jco",
                                  addaEllipses = TRUE,
                                   col.var = "contrib",
                                  gradient.cols = "BuGn",
                                  legend.title = list(fill = "Sexe", color = "Contrib"))

## Warning: ggrepel: 39 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

PCA issu du Package mixOmics

Pourquoi le graph plotInd de mixOmics semble être symétrique par rapport à celui de factominer L’individu est à-20 sur l’axe1 alors qu’il est à +20 avec ACP de factominer

trans.pca <-  pca(log10_CPMG, ncomp=5, center = TRUE, scale = TRUE)

plotIndiv(trans.pca, comp = c(1,2),ind.names = FALSE)

#CPMG_filt <-  CPMG[-match(c("L31R124"), table = rownames(CPMG)), ]


#PCA(CPMG, scale.unit = TRUE, ncp = 5, graph = TRUE)

Sparse PCA

Regarde si on peut observer naturellement une clusterisation en fonction d’une catégorie ? Quelque soit la catégorie sexe, age tabac, nous n’observons aucune clusterisation

res.spca.CPMG <- spca(log10_CPMG, scale=TRUE,
                         ncomp=3, keepX=c(10,10,10))

plotIndiv(res.spca.CPMG, ind.names = FALSE,
          group=Sample_ID$SEX,
          pch = as.numeric(Sample_ID$SEX),
          pch.levels =Sample_ID$SEX,
          legend = TRUE)

PLS-DA

Pourquoi l’axe X variate1 est plus faible que l’axe2. Possibilité de changer manuellement avec comp=c(2,1) X.label

trans.plsda <-  plsda(log10_CPMG, Sample_ID$SEX ,ncomp=5, scale = TRUE)

plotIndiv(trans.plsda, ind.names = FALSE, legend=TRUE,
          ellipse = TRUE, title="PLS-DA sur CPMG")

Cercle de correlation PLS-DA

N’est représenté uniquement les buckets qui contribuent à au moins 70%

plotVar(trans.plsda, var.names = FALSE, cutoff = 0.7,
        title = "Correlation circle plot avec cutoff 0.7")

Test de prédiction

Possibilité de faire de la prédiction dans le cas ou les ind seraient en sous groupes, ce qui n’est pas notre cas. ceci est juste un test

trans.testplsda <-  plsda(log10_CPMG, Sample_ID$SEX ,ncomp=5, scale = TRUE)

background <- background.predict(trans.testplsda, comp.predicted=2,
                                dist = "max.dist") 

plotIndiv(trans.testplsda, ind.names = FALSE, legend=TRUE,
          ellipse = TRUE, title="PLS-DA sur CPMG", background = background)

Contribution des variables

Les variables sélectionnées pour la composante 1 sont exprimées de manière positives pour le sexe1 et négative pour le sexe2

MyResult.plsda2 <- plsda(log10_CPMG, Sample_ID$SEX, ncomp=3)

plotLoadings(MyResult.plsda2, contrib = 'max', method = 'mean')

Clusturisation par PLS-DA

#Ceci est un test pour faire de l'intégration. X = métabo et Y =OTU.
X <- log10_CPMG
Y <- log10(NOESY + 1e-8)   ### Il faudrait mettre OTU en Y pour commencer intégration.
MyResult.spls <- spls(X,Y, keepX = c(25, 25), keepY = c(25,25))  
plotIndiv(MyResult.spls)

png(file = "out.png", width = 800, height = 700)

fig <- cim(MyResult.spls, comp = 1)

dev.off()

## png 
##   2

Création du Dataset mixomics

ncomp correspond: The number of components to include in the model for each block (does not necessarily need to take the same value for each block). Je ne sais pas. J’ai mis 3 pour le tabac Ne fonctionne pas pour CPMG

Data_mix <- mixOmics(X=NOESY, Y=Sample_ID$tabac, ncomp=5)

## a Partial Least Squares - Discriminant Analysis is being performed (PLS-DA)

Data_mix1 <- mixOmics(X=CPMG, Y=Sample_ID$tabac, ncomp=5)

## a Partial Least Squares - Discriminant Analysis is being performed (PLS-DA)

#les 2 ensembles ne fonctionnent pas
#Data_mix_tot <- mixOmics(X=list(Noesy,Cpmg), Y=Sample_ID$tabac, ncomp=3)
#Message : 
# Error in Check.entry.wrapper.mint.block(X = X, Y = Y, indY = indY, ncomp = ncomp,  : 
#   Each block of 'X' must have a unique name.


#head(Sample_ID$tabac)
#help(package='mixOmics')

Données NOESY

• Non inclus car donne les mêmes choses que CPMG